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We solve numerically the problem of finding the potential and electric 
field around a negatively charged metallic contact on the surface of an 
n-type semiconductor. The semiconductor, which has permittivity c, , fills 
the half-space y < 0. The contact is an infinitely long strip of width 
2a, defined by y = 0, ^ x ^ 2a, — w <z< «. The region y > is 
vacuum with permittivity e . In suitable dimensionless coordinates the 
potential <f> satisfies Laplace's equation in y > and the equation V <f> = 
e* — 1 iny < 0. On the boundary y = 0, <f> = <f> < 0, ^ x ^ 2a, and 
the usual electromagnetic boundary conditions at the remainder of the 
interface. Finite difference schemes are used to solve the resulting boundary 
value problem. 

In most practical cases |<£ | » 1 and tj = e /ei <3C 1. We examine in 
considerable detail the limiting case tj = 0, first for the less practical 
situation where \<f> \ « 1 and then for |</» | » 1. In case the |«fr>| « 1 we show 
that our numerical solution agrees well with the exact analytical solution 
of a linearized version of the problem. For \<f> \ y> 1, we give plots of the 
equipotential curves, curves of equal charge density, and curves of constant 
electric field amplitude. These results also yield expressions for the capaci- 
tance of both a strip and a circular electrode. The modifications of these 
results when rj > are also given in some detail. Finally, we discuss the 
numerical calculations at some length. 

I. INTRODUCTION AND FORMULATION OF THE PROBLEM 

In the study of a number of solid-state devices, it is important to 
know the electrostatic potential in the neighborhood of a metal-semi- 
conductor contact (Schottky diode) and in a metal insulator-semi- 
conductor structure (MOS capacitor). 

Motivated by this interest, we consider in this paper the following 

f On leave from the Technion-Isreal Institute of Technology, Haifa, Israel, when 
this work was conducted. 
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problem. An infinitely long metallic strip of width 2a* and zero thick- 
ness occupies the region ^ x* ^ 2a*, y* = 0, as shown in Fig. 1. 
The region y* > is filled with air and the region y* < is filled with 
an rz-type semiconductor. The metallic strip is charged to a negative 
potential, <p% < 0, and we wish to calculate the potential, the electric 
field and the electric charge density distribution in the semiconductor 
under the assumption of no current flow. As will be shown later, this 
solution also determines the potential, field and charge density around 
a circular metallic contact. 

For the doping levels normally encountered in such devices, say 
10 20 -10 22 ra~ 3 the electrostatic potential <j>* in the semiconductor is 
governed by the Poisson equation, 



'*2,* _ 



= -p*A. , 



(i) 



where V* 2 is the Laplacian, e : is the permittivity of the semiconductor, 
and the net volume charge density p* is given by 1 



p* = qN d [l ~ exp (q<t>*/kT)\, 



(2) 



where — q is the charge of an electron, N d is the donor number density, 
k is Boltzmann's constant, and T is the absolute temperature. In 
writing down equation (2) we have neglected the contact potential ^ 
between the metal strip and the semiconductor, 2 since in many cases 
of practical interest | <p | ^> 1*1- Here and in the following, starred 
quantities have rationalized MKS dimensions; unstarred quantities, 
except for a few obvious physical parameters, are dimensionless. 

As in Ref. (1) we introduce the dimensionless length and potential 
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Fig. 1 — A diagram of the geometry of the problem PI showing the conducting 
strip at the air-semiconductor interface and the coordinate system used. The dimen- 
sionless parameters shown are defined in equation (3). The symbol [ ] is defined 
in equation (8). 
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by the relations 

x = x*/\ D , y = y*/\ D . « - «*A/> . (3) 

= q<f>*/kT, p = p*/qN d , 

where the Debye length is given by 3 

\d = (e l kT/ q 2 N,) i . (4) 

Typical values for a lightly doped semiconductor device are N d = 
10 21 m~ 3 , eo/e, = .0625 for germanium and e /ei = .078 for silicon, 
and a* = 10" 4 m, where e is the permittivity of free space. Then, for 
example, for silicon at T = 300°K, \ D = 1.35 X 10' 7 m, q/kT = 38.8 
volt" 1 , and a = 740. 

In terms of these dimensionless quantities the boundary value 
problem to be solved for the potential can be summarized mathe- 
matically as follows: 

(i) In the air, y > (see Fig. 1), the potential satisfies Laplace's 

equation 

V 2 = 0. (5) 

(it) In the semiconductor, y < 0, the potential satisfies Poisson- 
Boltzman's equation 

V 1 * - e* - 1. (6) 

(m) On the plate, ?/ = and ^ a; ^ 2a, the potential is a given 
(negative) constant 

= 0o < 0. (7) 

(iv) At the interface, y = 0, \x - a\ > a, the potential and the 
normal component of the electric displacement vector are continuous' 

[0(z, 0)] - 0(x-, 0+) - 0(a-, 0-) = 0, (Sa) 

[^H-'^™-^- -^ - (sb) 

where 

J] = Co/*. . ( 8c ) 

(v) At infinity 
where r = (a;" + y 2 ) h - 



lim (0) = 0, (9) 



r-*oo 
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We refer to the boundary value problem defined by conditions 
(i)-(v) as PI. We have been unable to solve this problem analytically, 
and instead have studied its solution numerically. 

Qualitatively, for large | 4> I the potential can be described easily. 
Within the semiconductor, the negative charge on the conducting 
strip repels the mobile electrons. This action produces a layer around 
the strip, called the depletion layer, from which almost all the mobile 
electrons have been expelled. The positive donor ions left behind make 
this a region of uniform positive volume charge density. Far from the 
plate the semiconductor is almost neutral, and these two regions are 
connected by a transition layer. When 77 = this transition layer is 
sharp and well defined and is several Debye lengths thick. However, 
when 77 5^ 0, this transition region becomes broader and diffuse near 
the semiconductor-air interface. In the air on the other hand, the 
potential is essentially that due to the dipole formed by the negative 
charge on the conducting strip and the positive charge due to the donor 
ions in the depletion layer. 

In Section II, we consider the special case y = 0. For bias voltages 
</>o which are typically encountered in semiconductor devices, the thick- 
ness of the depletion layer is large compared to a Debye length but 
small compared to the width of the strip, 2a. We will show that we 
then only need to consider a strip completely embedded in an n-type 
semiconductor. This will reduce the solution of PI to the solution of P2 

(i) In the semiconductor V 2 = e* — 1. (6) 

(it) On the plate, y = 0, ^ x ^ 2a, <t> - <f> < 0. (10) 

iiii) At infinity linw fa) = 0. (11) 

Solutions of P2 are obtained by the method of finite differences. For 
| <j) | « 1, equation (6) can be linearized to 

V 2 <£ = <f>. (12) 

The linearized version of P2 in the limit a = «> , has been solved ana- 
lytically by Lewis 6 and for small <f> his results agree excellently with 
our numerical solution in the region x < a. This provides a good check 
on our numerical methods. Note, however, that in our formulation we 
neglected the contact potential and therefore the problem is not physi- 
cally meaningful for this limit. For | <f> | » 1 numerical calculations of 
both the electric field and the potential are presented in considerable 
detail. Finally, the capacitance per unit length of the strip and the 
capacitance of a circular electrode are presented with particular em- 
phasis on edge effects. 
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In Section III we briefly discuss what modifications of the results of 
Section II must be made when y t± 0. It is shown that the chief effect 
of positive y in the semiconductor is a smearing out of the transition 
region near the air-semiconductor interface. 

In Section IV we give some of the details of the methods of numerical 
analysis used. 

II. THE CASE 77 = 

In most practical applications e„ « e 1 , as was pointed out in 
Section I, and it is therefore of interest to treat first the simpler, limit- 
ing problem with 17 = 0. When 77 = condition (8b) becomes 

t*(*,0-)=0, (13) 

and the solution for the potential in the semiconductor is decoupled 
from the solution for the potential in the air. In fact, the solution of 
PI in the semiconductor is now identical with the solution of P2 be- 
cause of symmetry. 

There are three characteristic lengths in the problem of the finite 
width strip: the half width of the strip, the Debye length, and the 
thickness of the depletion layer. If the half width of the strip is very 
large compared to both the Debye length and the thickness of the 
depletion layer, the solution below the plate and sufficiently far from 
the edges must be the one dimensional solution (independent of x). 
Thus, for a "sufficiently wide" finite strip, there are really only two 
characteristic lengths for the solution in the semiconductor: the Debye 
length and the depletion layer thickness. 

We determine more precisely what "sufficiently wide" means. It is 
well known that if one considers the one dimensional problem of an 
71-type semiconductor filling the region y < 0, with the plane y = 
held at a large negative potential <f> < 0, then the thickness of the 
depletion layer, R, is accurately given by 2 

R = I 20 o I*. (14) 

[The accuracy of this approximation is also discussed in Ref. (1).] The 
Debye length is equal to unity in our nondimensional coordinates. 
Thus, if in addition to -q = 0, we have 

a»l, (15a) 

and 

a»#, (15b) 
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then a disappears from the problem in the semiconductor as a char- 
acteristic length. Reverse bias voltages of the order of <f>$ ~ —5 to 
— 50 volts, (0o ~ —200 to —2000) are typical, and this corresponds to 
R ~ 19.7 to 63.0 at room temperature. Thus in most practical situa- 
tions condition (15) is satisfied, since a > 740, and we can expect the 
solutions to be independent of a also near the edges of the strip. 
If in addition to equations (13) and (15a) we have 

I 0o | « 1, (16) 

then P2 can be linearized, for equation (5) can be approximated by 

V 2 = 0, (17) 

since | | ^ | O | everywhere. This linear problem with a = <x> has 
been solved analytically by Lewis, 5 and in polar coordinates (see 
Fig. 2) the solution is given by 

0/0o = \ exp (r sin 6)\ 1 — erf r*( cos - -f sin -J \f 



+ \ exp (— r sin 9)\ 1 -f- erf r 



cos- -sin - 



(IS) 



where erf (z) is the error function. We have solved the (nonlinear) 
P2 by a finite difference method to be described in Section IV. The 
numerically calculated equipotentials are shown in Fig. 3, and in 
Fig. 4 we compare the numerical solution of P2 (crosses) with the 
analytical solution (18) of the linearized problem (continuous line). 



SEMICONDUCTOR ( t 





=e - 




• o < 



SEMICONDUCTOR ((, 

v 2 =e*-i 



Fig. 2 — A diagram of the geometry of the problem P2 showing the coordinate 
system used. 
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Fig. 3 — Equipotential curves for the case v = 0, <£o = —.01, and a = 4. 

These figures correspond to 4> Q = -.01 and a = 4, using the mesh 
size Ax = Ay = .1. The three curves in Fig. 4 are the potential along 
the lines b - b(x ^ 0, y = 0), c - c(x = 0, y ^ 0) and d - d(x = 4, 
7/ <; 0). These three lines are shown in Fig. 2. It is seen that the two 
solutions agree very well. 

In most practical cases, however, <f> is very large (typically 200 < 
| 4> 1 < 2000), so we shall concentrate on the large potential problem 
from now on. When tf> is not small, the problem cannot be linearized. 
We have been unable to find approximate analytic solutions, and so 
we have had to solve the problem numerically. 

From the results of Ref. (1), we should expect that if distances are 
normalized with respect to the depletion layer thickness, then the 
potential when normalized with respect to the plate potential should be 
essentially independent of <£ as | 4> \ —* °° • We introduce this normali- 
zation here: 

* = x/R, y = y/R, a = a/R, 4> = 4>/4>» , p = P, (19) 

where R is the depletion layer thickness given in equation (14). In 
terms of our new variables, the basic equation (6) for <$>(£ , y) reads 



V 2 = 



5? 



+ 



off 



= 2(1 - e* of ) = 2p. 



(20) 



We have solved P2 in the normalized (tilde) variables by the method 
of finite differences. A detailed description of the method will be given 
in Section IV. Calculations for <\> Q = —100 and —500 using the mesh 
sizes A£ = Ay = .05 and for <j> = — 2500 using the mesh sizes AX = 
Ay = .025 for a = 3 and 4, have been carried out. In Fig. 5 we show the 
equipotential (0) lines and in Fig. 6 the lines of constant charge p( = p) 
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Fig. 4 — A comparison of the analytical solution of the linearized problem 
(continuous line) with the numerical solution of the nonlinear problem P2 (crosses) 
for the case i\ = 0, # = —.01. The three curves are the potential along 6 — 6 
(* g 0, y » 0), c - c (a; = 0, y g 0) and d - d (x - 4, y ^ 0). 

for the case <f> = — 500. These curves show clearly the depletion layer 
and the transition layer. The curves for <6 = — 100 and — 2500 (which 
we do not show) are essentially the same; they differ only in that for 
<f> = — 100 the transition layer is thicker while for <j> = — 2500 it is 
sharper. It might be pointed out that in the tilde coordinates the 
thickness of the transition layer becomes vanishingly small as | <£ | — > ■» , 
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Fig. 5 — Equipotential curves for the case j? = 0, <j> a = —500, a = 3. 
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Fig. 6 — Constant charge density curves for the case tj = 0, <f> = -500 and a = 3. 



but in the true dimensional coordinates, the thickness is ■~4X £) , inde- 
pendent of 0o as | </>o | —> °° • 

Figures 5 and 6 also show that for £ > 2, the equipotential and equi- 
charge curves are essentially parallel to the £ axis. Thus a > 2R is a 
sufficient condition for replacing the strip of finite width by a semi- 
infinite strip. Furthermore, since d<p/d£ = d 2 $/d£ 2 = in the region 
£ > 2, the solution is one dimensional here. In Ref. (1) it was shown that 
an excellent approximate solution of the one dimensional problem is 
the "zeroth-order matching" solution 



= 



(1-11/ 

0, 



-l ^ y ^ 
y ^ -l. 



(21) 



From equation (21) we also see that for £ ^ 2 



E; = 

E s = 



d£ 

_d0 

dy 



(22) 



-2(1 - |£ 

o, 



-l £ £ 0, 

y ^ -l. 



In Fig. 7 we plot along the lines b - b(£ ^ 0, y = 0) c — c(£ = 
0, § ^ 0) and d — d(£ = 3, y ^ 0) for the case O = —500 and a = 3. 
We superimpose on this a plot of 0(3, y) as given by (21). The agreement 
on d — d between calculated by the finite difference method and 
given by (21) is excellent. This will be the case even for larger mesh 
sizes since the truncation error due to approximating the Laplacian by 
a five point difference scheme is proportional to the fourth derivatives 
of 0, which should be small since is essentially parabolic along d — d. 
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Fig. 7— Graphs of <£ along the lines b - b (x £ 0, ft = 0), c - c (£ = 0, $ ^ 0) 
and d — d (x = 3, y ^ 0) for the case 77 = 0, <£ = —500, and 5 = 3. Superimposed 
on the d — d curve is a plot of <^(3, y) (circles) given by equation (21). 

In Fig. 8 we plot | E s | along b — b and | E s | along c — c and d — d 
for the same case, and we superimpose on this a plot of | E e (3, y) \ 
as given by equation (22). 

Because of the singularities of the electric field (and all higher deriva- 
tives of the potential) near the edge £ = # = 0, one cannot expect to 
obtain a uniformly valid numerical solution there. Nevertheless, com- 
parison of the numerical and analytical solutions for the small potential 
case (Fig. 4) shows that even at the nearest mesh points to the plate 
edge, the error in the numerical computation of the potential is not 
large. The error in calculating the electric field is naturally greater. 
In order to decrease the truncation error (the difference between the 
exact solution of the differential equation and the solution of the dif- 
ference equations), one can decrease the mesh size uniformly over the 
computational field. This can become quite expensive and is not neces- 
sary. A more efficient scheme is to refine the mesh size only in a small 
region around the plate edge and to find the numerical solution inside 
this region using at the boundary the values obtained from the coarser 
grid. We have done this for a region of size 1X1 around the plate edge 
with a mesh size of A£ = A§ = .0125 in the case O = —500 and 5 = 3. 
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Fig. 8— A plot of \S g \ along b - b(£ g 0, y = 0) and \S S \ along c - c(£ = 0. 
7/ ^ 0) and d - d (£ = 3, ?/ ^ 0) for the case i» = JL 0o = -5(10, and fl = 3. 
Superimposed on the d-d curve is a plot (circles) of \E S (Z, y)\ given by equation 
(22). 

For this refined solution we plot along & — 6 and c — c in Fig. 9, and 
in Fig. 10 we plot | E £ | along b — b and | E Q \ along c — c. The squares 
are points obtained from the solution using the coarser mesh (A£ = 
by = .05), and we see that several grid points away from the plates, 
the two solutions agree nicely. 

In addition to equipotential curves, curves of constant field amplitude 
can also be plotted. We define 



[« 



ft*. s> = 5 + ss 



30 
5i// _ 



= HE 



(23) 



From equation (23) we see that max„- #(«, y) = $(£, 0) = 2 for * > 2. 
Near the plate edge, however, there are much higher fields. In Fig. 11 
we draw the curves of constant i£, again for the case <j> = —500 and 
a = 3. The plate edge region can be defined as the region where $(£, y) 
^ 2. In Fig. 12 we show in more detail the curves of constant $ in the 
plate edge region. The data is taken from the calculation of the potential 
using a refined mesh around the plate edge discussed above. 

An important application of the numerical solution is the computation 
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Fig. 9— Plots of along b - b (x ^ 0, y = 0) and c - c (£ = 0, y g 0) for the 
case tj = 0, 0o = —500, and a = 3. The solutions were obtained with a mesh size of 
A£ = Ay = .0125 (continuous line) and Ax = Ay = .05 (squares). 



of the capacitance per unit length of the strip, C*, where 



C* = 



dQ* 



(24) 



and Q* is the total (dimensional) charge per unit length on the strip. 
To calculate Q* we note that it is just equal and opposite to the total 
net charge per unit length in the semiconductor: 



Q* = - f f p*( x *, y*) dx* dy*. 



(25) 



We introduce the nondimensional (tilde) coordinates normalized to the 
depletion layer thickness and write 

Q* = -qNXoR 2 f f° pd£dy = 2e^ f f" pdsdy. (26) 

•'-oo •'-00 J— oo •'—00 

This last integral, which is almost independent of <j>* for large \<f>* \, 
was computed numerically by evaluating the sum ^ ^ p A£ Ay and 
we write the result in the form 
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Fig. 10— Plots of \E t \ along b - b(£ < 0, y = 0) and \Sg\ along c - c(£ = 0, 
# g 0) for the case jj = 0, 4> = —500, and a = 3. The solutions were obtained with 
a mesh size of A£ = Ay = .0125 (continuous line) and AX = Ay = .05 (squares). 
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Fig. 11— Curves of constant field amplitude for the case v = 0» *o = -500, 
and a = 3. 
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Fig. 12 — A detailed plot of the constant field curves near the edge of the plate 
for the case t) = 0, <j>o = —500 and a = 3. 



// 



P dx dy = 2a + 2D . 



(27) 



The first term on the right of equation (27) is the value we would obtain 
if the effects of the edges and the transition layer were ignored, 
and Z) is the correction for these effects. For a ^> 1 and | <f> \ ^> 1, D 
should be essentially independent of a and <f> , that is, of a* and <f>$ . 
For our basic computation, d = 3, O = —500 and Ax" = Ay = .05, we 
obtained 



D = 0.354. 



(28) 



This value of D was very insensitive to increasing d to 4, and it changed 
only slightly when | O I was increased to 2500 (with A£ = A-g = .025). 
When we insert equation (27) into (26), expressing a in dimensional 
coordinates, we obtain 



where 



Q* = 46,0*^— + D () ) , 
R* = Xo # = (-2 £ ^*/qN d ) h 



(29) 



(30) 



is the dimensional depletion layer thickness. From equations (24), (29) 
and (30) we get 



<"-^(*+»*S- 



(3D 
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It is clear that our calculations also give an accurate approximation to 
the potential around a conducting disc of radius r£ , situated at the air- 
semiconductor interface, and charged to a large negative potential, as 
long as 

r„ » R, (32) 



where 

r - r*/\ D . (33) 

In all previous expressions for the potential and fields we need only 
interpret x* as the radial coordinate r*. 

We can now also calculate the capacitance of a circular disc. In this 
case we have 



and 



where 



Q* = -qN d \lR* [ f f pdXdydz, (34) 

J— go •' — oo *' — 00 

f r /"* p d£ dij dz = irrl + 2irr D , (35) 

J— 00 J — 00 ■»— 00 



h = r /B, (36) 

and of course we still have D = 0.354. If we substitute equations (35) 
and (36) into (34) and differentiate with respect to <f>* , we get 



-R^{ 1 + 4A) r* 



C»-SS- l + 4D.«r • (37) 



Finally, the potential in the air, y > is related to the potential on 
the interface $(£, 0) by Green's formula 7 

4>(£, y) = -T yK* - & + fTT%. 0) #, (38) 

7T «/ — oo 

where implicit use has been made of boundary condition (8a). Using the 
values of <£(£, 0) obtained from the finite difference solutions in y < 0, 
the integral in (3S) has been evaluated numerically to give the potential 
in y > 0. 

III. THE CASE 77 7^ 

When 77 = eoAi is not zero, the problem of solving PI is complicated 
by the fact that the solutions for the potential in the air and in the 
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semiconductor become coupled and cannot be found separately. Never- 
theless, the effects of this coupling can be taken into account in the 
semiconductor, without actually solving the problem in the air, by 
using Green's formula (38). We did this and found the solution in the 
semiconductor by an iterative finite difference scheme. Using this solu- 
tion on the boundary, the integral in (38) was evaluated numerically 
to obtain the solution in the air. The details of this numerical scheme 
are given in Section IV. We continue to concentrate on the case of large 
potential (\<f> \ >>> 1) and a wide conductiong strip (a 2>> R), so the 
normalized (tilde) variables introduced in equation (19) are retained. 

Except in a wedge shaped region at the interface at each edge of the 
strip (see Fig. 13), the solution in the semiconductor is very insensitive 
to changes in 77 and is little different from the solution for rj = 0. The 
main new feature of the solution in the semiconductor is the appearance 
of a "shoulder" in the depletion layer near the surface, and a smearing 
out of the transition layer there. In fact, the larger 77 is, the less sharp 
the transition layer in this region is. This effect increases the charge in 
the wedges and thus also increases the capacitance of the plate. The 
solution in the air is, however, more sensitive to changes in 77. 

We illustrate these qualitative remarks by a number of graphs. The 
case illustrated in all these graphs is <f> = —500, a = 3, and the mesh 
size used was A.f = Ay = .05. In Fig. 14 we show the equipotential 
(<£) curves and in Fig. 13 we show the curves of constant charge (p) 
for 77 = .1. These two figures illustrate the previous remarks. In Fig. 15 
we give graphs of <£ along the lines b — b(£ S 0, y = 0) and c — c(£ = 
0, y ^ 0) for 77 = 0, .05, and .1. In Fig. 16 we give a vector plot of the 
field around the plate for 77 = .1. 

Finally, we have numerically evaluated the integral (27) in order to 
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Fig. 13 — Constant charge density curves for the case 77 = .1, O = —500, and 
a = 3. 
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Fig. 14 — Equipotential curves for the case rj = .1, <£ = — 5C0, and 3 = 3. 

calculate the capacitance for 77 ^ 0. This was done for O = —500, 
a = 3, for 77 = 0, .05, and .1 with AX = Ay = .05, and for 77 = 0, .05, 
.1, .15 and .2 with A£ = Ay = .1. For this range of values of 77, it was 
found that (27) and the expressions (31) and (37) for the capacitance 
remain valid if D = .354 is replaced by 



D(r,) = .354 + .43 i). 



(39) 



We believe the number .43 in equation (39) is correct to about 10 percent. 
If we had approximated the depiction layer by a rectangle with sides 
2a* and R* completed at each end with a quarter circle of radius R*, 
we would have obtained D(t?) = ir/4. With this value of D(y), equation 
(37) yields the approximation of Goodman 8 and of Sze and Gibbons. 
By a completely different technique Copeland 10 has independently 
obtained the values £(.078) = .375 and £(.0625) = .365 which agree 
rather well with our results. 



IV. DETAILS OF THE NUMERICAL METHODS 

In this section we discuss some of the details of the various finite 
difference schemes used in solving numerically the boundary value 
problems PI and P2. In each of the three cases studied, tj = and 
I 0o I « lj and 77 ^ 0, I 0o I » 1, the basic nonlinear partial differential 
equation (6) or (20) is replaced by a system of finite difference equations, 
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Fig. 15 — Graphs of along the lines b — b (x g 0, y = 0) and c — c (£ = 0, 
y ^ 0) for the case <£ = —500 and a = 3 for y = 0, .05, and .1. 

but the boundary conditions, or the method of solving the finite dif- 
ference equations differ from case to case. 

In each case the infinite plane is replaced by a finite rectangle, as 
shown in Fig. 17, where advantage is taken of the symmetry of the 
problem around x = a. The rectangle is subdivided into a square mesh 
with mesh spacing h, so the rectangle has length Mh and depth Nh. 
In the small potential case, we define 



4n ti = <t>(a — (i - l)/i, —(j - l)h, 
while in the large potential case 

*«,, = 4>{a - (i - l)/i, -(j - l)h), 



(40a) 



(40b) 



where (1 ^ * £ M + 1, 1 ^ i ^ iV + 1). In either case, at each interior 
mesh point, V 2 ^,,, is approximated by the five point formula 

Vfo., = (* +1 ., + *,--,., + 0,-. J + i + *«.,-i - 40,,,)A 2 . (41) 

At the boundary points (i = 1, 2 ^ ; ^ N) we make use of the basic 
symmetry to write 



Vfo,.,- = {2fa ti + # M+1 + $,,,-_, - 44>i,,)/h a , 



(42) 
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and similarly in the case 77 = we write 

VJfo,, = (0,- + ,.i + 0,-,,, + 20,. a - 40..,)A a (43) 

for (M, + 1 < i ^ M), where il^/i = a for | O | « 1 and MJi = a 
for |0 O I » 1. 

When 77 = 0, I 0o I « li equation (6) is replaced by the difference 
equations 

Vfc w -«p(-|*„|)-l ( 44 ) 

for (1 ^ i ^ M, 2 <: j ^ AT) and (Af a + 1 < i ^ M, ; = 1). Since 
^ everywhere, the replacement of exp by exp — | <f> | changes 
nothing analytically, but numerically it eliminates certain instabilities 
in some iteration schemes. Equation (43) must be supplemented by 
further equations obtained from the boundary conditions. The con- 
dition = 0o on the strip yields 

0,.,, =0 O , (1 = i = Af, + 1). (45) 
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Fig. 16 — A vector plot of the field around the edge of the plate for the case »? = .1, 
O = —500, and a = 3. Each line segment has the direction of the field and the 
length of the segment is proportional to the magnitude of the field. 
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Fig. 17 — A schematic drawing of the square grid used in the numerical calculations. 

On the remaining boundaries condition (11) is approximated by re- 
quiring the vanishing of the normal derivatives, which decay faster 
than the potential itself. This yields the finite difference equations 



i.i — 4>* 



(1 ^ 3 ^ N), 
(1 ^ % ^ M). 



(46) 
(47) 



Equations (44), (45), (46) and (47) form a system of MN + M + N 
equations for the values of </>,,, at MN + M + N mesh points fou+i.ir-H 
does not appear in these equations). 

This set of transcendental equations has a unique solution 11 which 
was determined numerically by a point successive over-relaxation itera- 
tive method. 12,13 Briefly described, the method is as follows: an initial 
guess, 0j") , is made for <£,,,. At each interior point subsequent iterations 
0,-"} are determined by the equations 

x(n + l) _ l/.fn) . ,(n+l> I . 00 j_ 

<Pi.i — 4l9«+i,/ + 9.-1. / + 0.-.; + i + 



(n + l) 
>».)-l 



— h 2 exp (— | (f>i 

.(n+l) T(n + D I /i \,l») 

<Pi.i = «9.-.y + (1 — aWi.i • 

The over-relaxation parameter co is given by 13 
a> = 1 + a 2 /[l + (1 - a 2 )*? 
<r = i{cos (tt/AO + COS (vr/il/)}. 



+ /i 2 }, 



(48) 
(49) 

(50a) 
(50b) 



This value of a> is optimal for the Dirichlet problem for Laplace's 



ELECTROSTATIC POTENTIAL 873 

equation. Equation (48) must be modified in an obvious way at the 
boundary points where (44) is satisfied, and the <£,•,,• at the remaining 
boundary points must be eliminated with the aid of equations (45) 
through (47). The iterations are repeated until either 

5 M = [maxU?; 1 ' -fpiSD/l*! (51a) 



or 



5 2 =(EZI *57" " <*>" HVll *« i (MN + M + N)*\ (51b) 

is suitably small. It is important to note that the nonlinear term, 
exp (<f>i,i), appears on the right hand side of equation (48). For this 
reason the method is referred to as an explicit scheme. Typically in the 
small potential case we chose M, = 40, M = 100, N = 92, h = A, and 
0j.°) = 0, and after 50 iterations we obtained 5 W < 10" 3 and 8 2 < 10~ 4 . 
The iteration scheme just outlined is only conditionally stable, and 
it can be shown that an approximate condition for its convergence when 
applied to a system of equations of the form 

Vfa,,, = /(*..,) (52) 

is 

«(l + lj j) < 2, (53a) 

where 

C = max /'(<£). (53b) 

When | 0o | « 1) condition (53a) does not impose a severe restriction, 
since the region over which the solution varies significantly is small. 
Then h, M and N can be chosen so that M and N are not too large 
and at the same time h and w satisfy equation (52) for economically 
acceptable w. When | <£ | » 1, however, the depletion layer is so large 
that M and iV must be chosen excessively large in order for (52) to 
be satisfied. 

In order to avoid this difficulty, we have employed an implicit scheme 
instead of an explicit one. For convenience we first introduce rescaled 
(tilde) coordinates in (19). For y = 0, the resulting finite difference 
equations are just those given in equations (41) through (47) except that 
(44) is replaced by 

Vfr,-., = 211 - exp (-|0„0„ |)|. (54) 
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The explicit iteration scheme defined by equations (48) and (49) is 
replaced by the following implicit scheme. The quantity <?-", +1) is now 
the root of the transcendental equation 

$iT -f «p Marts 11 I) 

= *wa.i + *.--t!) + fR« + *, { :£i - 2^ a j . (55) 

The solution of this equation, which can be found by the Newton- 
Raphson method, is then substituted into equation (49) to obtain 
<£i"t u . Equation (55) must be modified in an obvious way at the bound- 
ary points where (43) is satisfied, and the fa , ,• at the remaining boundary 
points are eliminated as before. When < w < 2, the iteration scheme 
just defined can be shown to be unconditionally stable, 12 that is, it 
converges for any mesh size h and any 0j°) . However, in practice we 
found a) ~ 1.5 to provide the most rapid convergence. The rapidity of 
convergence is not very sensitive to small changes in a> around a> = 1.5. 
Typically in the large potential case for 77 = we chose M x = 60, 
M = 100, N = 40, and h = .05, and after ~80 iterations we obtained 
5*. < 10" 4 and 5 2 < 1.6 X 10~ 5 . 

When rj 7^ 0, this iteration scheme must be modified to take into ac- 
count the coupling of the potential in the semiconductor to the potential 
in the air. If we attempt to solve directly for the potential in both the 
air and in the semiconductor by the method of finite differences, we 
encounter difficulties. The potential in the air decays very slowly at 
infinity, so that if we use a reasonable number of mesh points, boundary 
conditions such as (46) and (47) introduce fairly large errors into the 
calculation. In order to circumvent this problem without using an 
inordinately large number of mesh points we proceed as follows. We 
replace boundary condition (8b) by the equivalent condition that 4 

~L'B dl * = LI >* ****** (56) 

where G is a small square shown in Fig. 17, of side h centered at the 
boundary point (M 1 + 1 < i < M, j = 1). The integral on the left of 
equation (56) is a line integral around the boundary of G, and d/dn* 
denotes the normal derivative. The finite difference approximation to 
(56) yields the equation in tilde coordinates 

</>,-, 2 + #,.0 + ^ * - (4>i + l,l + 0.-I.O - 2(1 + 7?)0,M 

= h 2 il - exp(-|<M»,-.i I)}, (57) 
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where <£,, is the potential at the mesh point in the air (y = h). Equation 
(57) now replaces (43) at the air-semiconductor interface points, and 
yields the iteration equations 

+ on Z S ll*S + *H " *'} - CM-, + K t < M). (58) 

Once the (n + 1) iterates 0<? f * n have been determined for all points 
within and on the boundary of the semiconductor, the^J^" , {M i + 1 < 
i fS M), are calculated from the finite difference version of Green's 
formula (38) using the values 0£J" already known. We did not prove 
that this modified iteration scheme converges, but it works well in 
practice. In this case also, the empirically determined value w = 1.5 
seems optimal. 

Finally, we make several brief comments about the difference between 
the true solutions of PI or P2 and the numerical solutions, that is, the 
truncation errors. It is well known 13 that when the true solution is 
smooth enough (specifically when the fourth derivatives are bounded) 
the truncation error is OQi 2 ). However, in our problem, near the plate's 
edge x = y = 0, even the first derivative is not bounded, and the above 
estimate fails. Wasow" considered a problem similar to ours with smooth 
boundaries and piecewise analytic boundary values and found that the 
truncation error vanished when h — > 0. 

It can be shown 15 that near the corner x = y = 0, the singularity in 
the field for arbitrary -q is 

V0 = 0(r~ h ) (59) 

for both the small and large potential cases (in the appropriate coordi- 
nates) where r is the distance from the corner. Bramble, and others, 16 
investigated the truncation (or discretization) error of such problems, 
and if we combine equation (59) with their Theorem 3.1, we find that 
near the corner the truncation error is 

I*,./ - *«..(*, y) I = 0(**), (60) 

and it vanishes as h —> 0. 

The truncation error is not uniform, and because of the slow con- 
vergence near the corner, it is worthwhile to refine the mesh size in a 
small region around it (as was described in Section II). This was done 
for the small potential case for tj = and (f> = —.01, for a refined 
mesh size of h = .0125 and in Fig. 18 we compare the analytic solution 
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Fig. 18 — A comparison of the analytical solution of the linearized problem (con- 
tinuous line) with the numerical solution of the nonlinear problem P2 (crosses) for 
the case ij = 0, <f>o = .01 and for a mesh size near the plate edge of Ax = Ay = .0125. 

(continuous line) with the numerical one (crosses) along the lines 6 — 6 
and c — c. 

All the calculations described in this paper were performed on a GE 
635 digital computer. The graphical output was obtained with the aid 
of a routine for calculating level curves written by G. S. Deem, L. K. 
Russell, and N. J. Zabusky. 17 
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